非線形回帰分析 - 17

指数関数のパラメータの推定

a0,をシフト

\(\Large \displaystyle y_i = a_0 \ Exp (- a_1 x_i) + a_2 \)

 

a0       9.1640 9.6640 9.9640 10.0640 10.1640 10.2640 10.3640 10.6640 11.1640
δ       -1.0 -0.5 -0.2 -0.1 0.0 0.1 0.2 0.5 1.0
a1       0.5321 0.5134 0.5021 0.4983 0.4946 0.4908 0.4870 0.4755 0.4559
a2       1.3123 1.1040 0.9735 0.9290 0.8839 0.8383 0.7921 0.6500 0.3993
                         
i x y   \( \hat{y} \)
1 0 10   10.4763 10.7680 10.9375 10.9930 11.0479 11.1023 11.1562 11.3140 11.5633
2 2 4   4.4739 4.5654 4.6236 4.6436 4.6639 4.6846 4.7055 4.7702 4.8850
3 3 2   3.1694 3.1756 3.1828 3.1858 3.1891 3.1928 3.1968 3.2111 3.2427
4 4 1   2.4031 2.3438 2.3107 2.3001 2.2897 2.2796 2.2698 2.2419 2.2017
5 6 0.5   1.6886 1.5481 1.4633 1.4350 1.4067 1.3784 1.3501 1.2651 1.1235
6 9 0.1   1.3886 1.1992 1.0821 1.0425 1.0025 0.9622 0.9216 0.7977 0.5837
         
S (\(y_i - \hat{y} \)の平方和)       0.8610 0.4039 0.2771 0.2591

0.2531
(Se)

0.2591 0.2769 0.4010 0.8384
dS (Seとの差分)       0.6078 0.1507 0.0240 0.0060 0 0.0060 0.0238 0.1479 0.5853
                         

 

・残差平方和

推定値からの残差

\(\Large \displaystyle Se = \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) - \hat{a_2} \right]^2 \)

a0をシフトさせたときの,推定値からの残差

\(\Large \displaystyle Se = \sum_{i=1}^{n} \left[ y_i -a_0 \ Exp(- \hat{a_1} \ x_i) -\hat{a_2} \right]^2 \)

であり,a0を,δ,だけシフトさせて,固定し,その際のa1, a2の推定値をソルバーで推定しました.

dS,を見ていただけるとわかるように,推定値,Seが一番小さく,ほぼ左右対称に増加していることがわかります.

グラフ化すると,

のように,二乗+定数できれいに近似できます.

二次曲線の近似においてもきれいに近似でき,

\(\Large \displaystyle y = 0.2532 + 0.5965 \ \delta^2 \)

ここで,分散値は,

・分散

\(\Large \displaystyle Ve = \frac{1}{n-3} \sum_{i=1}^{n} \left[ y_i -\hat{a_0} \ Exp(- \hat{a_1} \ x_i) \right]^2 = \frac{Se}{n-3} = \frac{0.2532}{6-3} = 0.08438 \)

であり(a0,a1,a2の二つのパラメータが3つあるので,自由度は,n-3),

\(\Large \displaystyle 0.5965 \ \delta^2 = 0.08438 \)

となるδがSEとなるので,

\(\Large \displaystyle \delta^2 = \frac{ 0.08438}{0.5965} =1.4145 \)

\(\Large \displaystyle SE_{a_0} = \sqrt{\delta^2} =\color{red}{0.3761} \)

と推定できます.

 

・Rによる推定

Rでの近似を行ってみると,

プログラムは,

xx <- c(0,2,3,4,6,9)
yy <- c(11,5,3,2,1.5,1.1)
plot(xx,yy)
fm<-nls(yy~a0*exp(-a1*xx)+a2,start=c(a0=10,a1=0.5,a2=1),trace=TRUE)
summary(fm)

で,結果は,

Parameters:
Estimate Std. Error t value Pr(>|t|)
a0 10.16402 0.37708 26.954 0.000112 ***
a1 0.49456 0.04408 11.219 0.001518 **
a2 0.88393 0.26931 3.282 0.046348 *

となり,Kyplotにおいても,

推定値 標準誤差(SE)
A1 10.16401 0.377081
A2 0.494562 0.044083
A3 0.883926 0.269308

と同じ結果となり,今回の推定値と,ほぼ一致,します.

微妙に異なるのが気になりますが....

「統計解析の初歩」,の「1.2 非線形最小2乗法の基本的な考え方」には,

δのずらした値0.5 を変えると結果は異なり、近似標準誤差の精度が変わる
統計パッケージにより、近似標準誤差の値は幾分異なる

とあります.ここで,”0.5”,がどこから出てきたかはわかりません.そもそも横軸(x軸)の範囲に依存しちゃいますし..

そこで,σをとる値(±での平均値)でどう推定値が変わるかを調べてみました.その結果が,

とどのδでもRなどの推定値より下回っていました....謎です...

次に,指数関数+baseの肩のパラメータ,a1について,確認してましょう.

 

l tr